{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare times between the test cases run using Python/C code (int2e-speed.ipynb) and Julia code. Here's the Julia code from pyquante.jl:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "gcf (generic function with 4 methods)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function gaussian_product_center(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,\n",
    "                                    bexpn::Float64,bx::Float64,by::Float64,bz::Float64)\n",
    "    return (aexpn*[ax,ay,az]+bexpn*[bx,by,bz])/(aexpn+bexpn)    \n",
    "end\n",
    "dist2(dx,dy,dz) = dx*dx+dy*dy+dz*dz # Is there something in the standard library that does this?\n",
    "function Fgamma(m::Int64,x::Float64,SMALL::Float64=1e-12)\n",
    "    #println(\"Fgamma($m,$x)\")\n",
    "    x = max(x,SMALL) # Evidently needs underflow protection\n",
    "    return 0.5*x^(-m-0.5)*gammainc(m+0.5,x)\n",
    "end\n",
    "\n",
    "function gammainc(a::Float64,x::Float64)\n",
    "    # This is the series version of gamma from pyquante. For reasons I don't get, it \n",
    "    # doesn't work around a=1. This works alright, but is only a stopgap solution\n",
    "    # until Julia gets an incomplete gamma function programmed\n",
    "    if abs(a-1) < 1e-3\n",
    "        println(\"Warning: gammainc_series is known to have problems for a ~ 1\")\n",
    "    end\n",
    "    if x < (a+1.0)\n",
    "        #Use the series representation\n",
    "        gam,gln = gser(a,x)\n",
    "    else \n",
    "        #Use continued fractions\n",
    "        gamc,gln = gcf(a,x)\n",
    "        gam = 1-gamc\n",
    "    end\n",
    "    return exp(gln)*gam\n",
    "end\n",
    "\n",
    "function gser(a::Float64,x::Float64,ITMAX::Int64=100,EPS::Float64=3e-9)\n",
    "    # Series representation of Gamma. NumRec sect 6.1.\n",
    "    gln=lgamma(a)\n",
    "    if x == 0\n",
    "        return 0,gln\n",
    "    end\n",
    "    ap = a\n",
    "    delt = s = 1/a\n",
    "    for i in 1:ITMAX\n",
    "        ap += 1\n",
    "        delt *= (x/ap)\n",
    "        s += delt\n",
    "        if abs(delt) < abs(s)*EPS\n",
    "            break\n",
    "        end\n",
    "    end\n",
    "    return s*exp(-x+a*log(x)-gln),gln\n",
    "end\n",
    "\n",
    "function gcf(a::Float64,x::Float64,ITMAX::Int64=200,EPS::Float64=3e-9,FPMIN::Float64=1e-30)\n",
    "    #Continued fraction representation of Gamma. NumRec sect 6.1\"\n",
    "    gln=lgamma(a)\n",
    "    b=x+1.-a\n",
    "    c=1./FPMIN\n",
    "    d=1./b\n",
    "    h=d\n",
    "    for i in 1:ITMAX\n",
    "        an=-i*(i-a)\n",
    "        b=b+2.\n",
    "        d=an*d+b\n",
    "        if abs(d) < FPMIN\n",
    "            d=FPMIN\n",
    "        end\n",
    "        c=b+an/c\n",
    "        if abs(c) < FPMIN\n",
    "            c=FPMIN\n",
    "        end\n",
    "        d=1./d\n",
    "        delt=d*c\n",
    "        h=h*delt\n",
    "        if abs(delt-1.) < EPS\n",
    "            break\n",
    "        end\n",
    "    end\n",
    "    gammcf = exp(-x+a*log(x)-gln)*h\n",
    "    return gammcf,gln\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "vrr (generic function with 1 method)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function vrr(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,aI::Int64,aJ::Int64,aK::Int64,\n",
    "        bexpn::Float64,bx::Float64,by::Float64,bz::Float64,\n",
    "        cexpn::Float64,cx::Float64,cy::Float64,cz::Float64,cI::Int64,cJ::Int64,cK::Int64,\n",
    "        dexpn::Float64,dx::Float64,dy::Float64,dz::Float64,m::Int64)\n",
    "    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)\n",
    "    qx,qy,qz = gaussian_product_center(cexpn,cx,cy,cz,dexpn,dx,dy,dz)\n",
    "    zeta,eta = aexpn+bexpn,cexpn+dexpn\n",
    "    wx,wy,wz = gaussian_product_center(zeta,px,py,pz,eta,qx,qy,qz)\n",
    "    #println(\"P: $px,$py,$pz, Q: $qx,$qy,$qz, W: $wx,$wy,$wz, $zeta,$eta\")\n",
    "    \n",
    "    val = 0\n",
    "    if cK>0\n",
    "        val = (qz-cz)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "            cexpn,cx,cy,cz,cI,cJ,cK-1,dexpn,dx,dy,dz,m) +\n",
    "            (wz-qz)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK-1,dexpn,dx,dy,dz,m+1)\n",
    "        #println(\"val1=$val\")\n",
    "        if cK>1\n",
    "            val += 0.5*(cK-1)/eta*(\n",
    "                vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "                    cexpn,cx,cy,cz,cI,cJ,cK-2,dexpn,dx,dy,dz,m) -\n",
    "            zeta/(zeta+eta)*\n",
    "                vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "                    cexpn,cx,cy,cz,cI,cJ,cK-2,dexpn,dx,dy,dz,m+1) )\n",
    "        #println(\"val2=$val\")\n",
    "        end\n",
    "        if aK>0\n",
    "            val += 0.5*aK/(zeta+eta)*\n",
    "                vrr(aexpn,ax,ay,az,aI,aJ,aK-1,bexpn,bx,by,bz,\n",
    "                    cexpn,cx,cy,cz,cI,cJ,cK-1,dexpn,dx,dy,dz,m+1)\n",
    "        end\n",
    "        #println(\"val3=$val\")\n",
    "        return val\n",
    "    elseif cJ>0\n",
    "        val = (qy-cy)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "            cexpn,cx,cy,cz,cI,cJ-1,cK,dexpn,dx,dy,dz,m) +\n",
    "        (wy-qy)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ-1,cK-1,dexpn,dx,dy,dz,m+1)\n",
    "        #println(\"val4=$val\")\n",
    "        if cJ>1\n",
    "            val += 0.5*(cJ-1)/eta*(\n",
    "            vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ-2,cK,dexpn,dx,dy,dz,m) -\n",
    "            zeta/(zeta+eta)*\n",
    "            vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ-2,cK,dexpn,dx,dy,dz,m+1)\n",
    "            )\n",
    "        #println(\"val5=$val\")\n",
    "        end\n",
    "        if aJ>0\n",
    "            val += 0.5*aJ/(zeta+eta)*\n",
    "            vrr(aexpn,ax,ay,az,aI,aJ-1,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ-1,cK,dexpn,dx,dy,dz,m+1)\n",
    "        end\n",
    "        #println(\"val6=$val\")\n",
    "        return val\n",
    "    elseif cI>0\n",
    "        val = (qx-cx)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "            cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m) +\n",
    "        (wx-qx)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI-1,cJ,cK-1,dexpn,dx,dy,dz,m+1)\n",
    "        #println(\"val7=$val\")\n",
    "        if cI>1\n",
    "            val += 0.5*(cI-1)/eta*(\n",
    "            vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI-2,cJ,cK,dexpn,dx,dy,dz,m) -\n",
    "            zeta/(zeta+eta)*\n",
    "            vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI-2,cJ,cK,dexpn,dx,dy,dz,m+1)\n",
    "            )\n",
    "        end\n",
    "        #println(\"val8=$val\")\n",
    "        if aI>0\n",
    "            val += 0.5*aI/(zeta+eta)*\n",
    "            vrr(aexpn,ax,ay,az,aI-1,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m+1)\n",
    "        end\n",
    "        #println(\"val9=$val\")\n",
    "        return val\n",
    "    elseif aK>0\n",
    "        val = (pz-az)*vrr(aexpn,ax,ay,az,aI,aJ,aK-1,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) +\n",
    "        (wz-pz)*vrr(aexpn,ax,ay,az,aI,aJ,aK-1,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)\n",
    "        #println(\"val10=$val\")\n",
    "        if aK>1\n",
    "            val += 0.5*(aK-1)/zeta*(\n",
    "            vrr(aexpn,ax,ay,az,aI,aJ,aK-2,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m) -\n",
    "            eta/(zeta+eta)*\n",
    "            vrr(aexpn,ax,ay,az,aI,aJ,aK-2,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m+1)\n",
    "            )\n",
    "        end\n",
    "        #println(\"val11=$val\")\n",
    "        return val\n",
    "    elseif aJ>0\n",
    "        val = (py-ay)*vrr(aexpn,ax,ay,az,aI,aJ-1,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m)+\n",
    "        (wy-py)*vrr(aexpn,ax,ay,az,aI,aJ-1,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)\n",
    "        #println(\"val12=$val\")\n",
    "        if aJ>1\n",
    "            val += 0.5*(aJ-1)/zeta*(\n",
    "            vrr(aexpn,ax,ay,az,aI,aJ-2,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) -\n",
    "            eta/(zeta+eta)*\n",
    "            vrr(aexpn,ax,ay,az,aI,aJ-2,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)\n",
    "            )\n",
    "        end\n",
    "        #println(\"val13=$val\")\n",
    "        return val\n",
    "    elseif aI>0\n",
    "        val = (px-ax)*vrr(aexpn,ax,ay,az,aI-1,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) +\n",
    "        (wx-px)*vrr(aexpn,ax,ay,az,aI-1,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)\n",
    "        #println(\"val14=$val\")\n",
    "        if aI>1\n",
    "            val += 0.5*(aI-1)/zeta*(\n",
    "            vrr(aexpn,ax,ay,az,aI-2,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) -\n",
    "            eta/(zeta+eta)*\n",
    "            vrr(aexpn,ax,ay,az,aI-2,aJ,aK,bexpn,bx,by,bz,\n",
    "                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)\n",
    "            )\n",
    "        end\n",
    "        #println(\"val15=$val\")\n",
    "        return val\n",
    "    end\n",
    "\n",
    "    rab2 = dist2(ax-bx,ay-by,az-bz)\n",
    "    rcd2 = dist2(cx-dx,cy-dy,cz-dz)\n",
    "    rpq2 = dist2(px-qx,py-qy,pz-qz)\n",
    "    T = zeta*eta/(zeta+eta)*rpq2\n",
    "    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)\n",
    "    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)\n",
    "    #println(\"rab2=$rab2,rcd2=$rcd2,rpq2=$rpq2,T=$T,Kab=$Kab,Kcd=$Kcd\")\n",
    "    return Kab*Kcd/sqrt(zeta+eta)*Fgamma(m,T)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "vrr_iter (generic function with 1 method)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function vrr_iter(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,aI::Int64,aJ::Int64,aK::Int64,\n",
    "        bexpn::Float64,bx::Float64,by::Float64,bz::Float64,\n",
    "        cexpn::Float64,cx::Float64,cy::Float64,cz::Float64,cI::Int64,cJ::Int64,cK::Int64,\n",
    "        dexpn::Float64,dx::Float64,dy::Float64,dz::Float64,M::Int64)\n",
    "    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)\n",
    "    qx,qy,qz = gaussian_product_center(cexpn,cx,cy,cz,dexpn,dx,dy,dz)\n",
    "    zeta,eta = aexpn+bexpn,cexpn+dexpn\n",
    "    wx,wy,wz = gaussian_product_center(zeta,px,py,pz,eta,qx,qy,qz)\n",
    "    rab2 = dist2(ax-bx,ay-by,az-bz)\n",
    "    rcd2 = dist2(cx-dx,cy-dy,cz-dz)\n",
    "    rpq2 = dist2(px-qx,py-qy,pz-qz)\n",
    "    T = zeta*eta/(zeta+eta)*rpq2\n",
    "    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)\n",
    "    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)\n",
    "    mtot = aI+aJ+aK+cI+cJ+cK+M\n",
    "\n",
    "    vrr_terms = zeros(Float64,(aI+1,aJ+1,aK+1,cI+1,cJ+1,cK+1,mtot+1))\n",
    "    \n",
    "    for m in 0:mtot\n",
    "        vrr_terms[1,1,1, 1,1,1, m+1] = Fgamma(m,T)*Kab*Kcd/sqrt(zeta+eta)\n",
    "    end\n",
    "    \n",
    "    for i in 0:(aI-1)\n",
    "        for m in 0:(mtot-i-1)\n",
    "            vrr_terms[i+2,1,1, 1,1,1, m+1] = (\n",
    "                 (px-ax)*vrr_terms[i+1,1,1, 1,1,1, m+1] + \n",
    "                 (wx-px)*vrr_terms[i+1,1,1, 1,1,1, m+2])\n",
    "\n",
    "            if i>0\n",
    "                vrr_terms[i+2,1,1, 1,1,1, m+1] += i/2/zeta*(\n",
    "                    vrr_terms[i,1,1, 1,1,1, m+1] -\n",
    "                    eta/(zeta+eta)*vrr_terms[i,1,1, 1,1,1, m+2])\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "#=\n",
    "    for j in 0:(aJ-1)\n",
    "        for i in 0:aI\n",
    "            for m in 0:(mtot-i-j-1)\n",
    "                println((\"b\",i,j,m))\n",
    "                vrr_terms[i+1,j+2,1, 1,1,1, m+1] = (\n",
    "                (py-ay)*vrr_terms[i+1,j+1,1, 1,1,1, m+1] +\n",
    "                (wy-py)*vrr_terms[i+1,j+1,1, 1,1,1, m+2])\n",
    "                if j>0\n",
    "                    vrr_terms[i+1,j+2,1, 1,1,1, m+1] += j/2/zeta*(\n",
    "                        vrr_terms[i+1,j,1, 1,1,1, m+1] -\n",
    "                        eta/(zeta+eta)*vrr_terms[i+1,j,1, 1,1,1, m+2])\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "\n",
    "    for k in 0:(aK-1)\n",
    "        for j in 0:aJ\n",
    "            for i in 0:aI\n",
    "                for m in 0:(mtot-i-j-k-1)\n",
    "                    println((\"c\",i,j,k,m))\n",
    "                    vrr_terms[i+1,j+1,k+2, 1,1,1, m+1] = (\n",
    "                    (pz-az)*vrr_terms[i+1,j+1,k+1, 1,1,1, m+1] +\n",
    "                    (wz-pz)*vrr_terms[i+1,j+1,k+1, 1,1,1, m+2])\n",
    "                    if k>0\n",
    "                        vrr_terms[i+1,j+1,k+2, 1,1,1, m+1] += k/2/zeta*(\n",
    "                        vrr_terms[i+1,j+1,k, 1,1,1, m+1] -\n",
    "                        eta/(zeta+eta)*vrr_terms[i+1,j+1,k, 1,1,1, m+2])\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "\n",
    "    for q in 0:(cI-1)\n",
    "        for k in 0:aK\n",
    "            for j in 0:aJ\n",
    "                for i in 0:aI\n",
    "                    for m in 0:(mtot-i-j-k-q-1)\n",
    "                        println((\"d\",i,j,k,q,m))\n",
    "                        vrr_terms[i+1,j+1,k+1, q+2,1,1, m+1] = (\n",
    "                        (qx-cx)*vrr_terms[i+1,j+1,k+1, q+1,1,1, m+1] +\n",
    "                        (wx-qx)*vrr_terms[i+1,j+1,k+1, q+1,1,1, m+2])\n",
    "                        if q>0\n",
    "                            vrr_terms[i+1,j+1,k+1, q+2,1,1, m+1] += q/2/eta*(\n",
    "                            vrr_terms[i+1,j+1,k+1, q,1,1, m+1] -\n",
    "                            eta/(zeta+eta)*vrr_terms[i+1,j+1,k+1, q,1,1, m+2])\n",
    "                        end\n",
    "                        if i>0\n",
    "                            vrr_terms[i+1,j+1,k+1, q+2,1,1, m+1] += (\n",
    "                            i/2/(zeta+eta)*vrr_terms[i,j+1,j+1, q+1,1,1, m+2])\n",
    "                        end\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "\n",
    "    for r in 0:(cJ-1)\n",
    "        for q in 0:cI\n",
    "            for k in 0:aK\n",
    "                for j in 0:aJ\n",
    "                    for i in 0:aI\n",
    "                        for m in 0:(mtot-i-j-k-q-r-1)\n",
    "                            println((\"e\",i,j,k,q,r,m))\n",
    "                            vrr_terms[i+1,j+1,k+1, q+1,r+2,1, m+1] = (\n",
    "                            (qy-cy)*vrr_terms[i+1,j+1,k+1, q+1,r+1,1, m+1] +\n",
    "                            (wy-qy)*vrr_terms[i+1,j+1,k+1, q+1,r+1,1, m+2])\n",
    "                            if r>0\n",
    "                                vrr_terms[i+1,j+1,k+1, q+1,r+2,1, m+1] += r/2/eta*(\n",
    "                                vrr_terms[i+1,j+1,k+1, q+1,r+1,1, m+1] -\n",
    "                                zeta/(zeta+eta)*vrr_terms[i+1,j+1,k+1, q+1,r,1, m+2])\n",
    "                            end\n",
    "                            if j>0\n",
    "                                vrr_terms[i+1,j+1,k+1, q+1,r+2,1, m+1] += (\n",
    "                                j/2/(zeta+eta)*vrr_terms[i+1,j,k+1, q+1,r+1,1, m+2])\n",
    "                            end\n",
    "                        end\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "\n",
    "    for s in 0:(cK-1)\n",
    "        for r in 0:cJ\n",
    "            for q in 0:cI\n",
    "                for k in 0:aK\n",
    "                    for j in 0:aJ\n",
    "                        for i in 0:aI\n",
    "                            for m in 0:(mtot-i-j-k-q-r-s-1)\n",
    "                                println((\"f\",i,j,k,q,r,s,m))\n",
    "                                vrr_terms[i+1,j+1,k+1, q+1,r+1,s+2, m+1] = (\n",
    "                                (qz-cz)*vrr_terms[i+1,j+1,k+1, q+1,r+1,s+1, m+1]+\n",
    "                                (wz-qz)*vrr_terms[i+1,j+1,k+1, q+1,r+1,s+1, m+2])\n",
    "                                if s>0\n",
    "                                    vrr_terms[i+1,j+1,k+1, q+1,r+1,s+2,m+1] += s/2/eta*(\n",
    "                                    vrr_terms[i+1,j+1,k+1,q+1,r+1,s, m+1] -\n",
    "                                    zeta/(zeta+eta)*vrr_terms[i+1,j+1,k+1,q+1,r+1,s,m+2])\n",
    "                                end\n",
    "                                if k>0\n",
    "                                    vrr_terms[i+1,j+1,k+1, q+1,r+1,s+2,m+1] += (\n",
    "                                    k/2/(zeta+eta)*vrr_terms[i+1,j+1,k, q+1,r+1,s+1,m+2])\n",
    "                                end\n",
    "                            end\n",
    "                        end\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    println(\"before return\")\n",
    "    =#\n",
    "    vrr_terms[aI+1,aJ+1,aK+1,cI+1,cJ+1,cK+1,M+1]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "test_vrr (generic function with 1 method)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function test_vrr()\n",
    "    ax=ay=az=bx=by=bz=cx=cy=cz=dx=dy=dz=0.0\n",
    "    aexpn=bexpn=cexpn=dexpn=1.0\n",
    "    aI=aJ=aK=0\n",
    "    cI=cJ=cK=0\n",
    "    M=0\n",
    "\n",
    "    for (ax,ay,az, aI,aJ,aK, cI,cJ,cK, result) in [\n",
    "            (0.,0.,0., 0,0,0, 0,0,0, 4.37335456733),\n",
    "            (0.,0.,0., 1,0,0, 1,0,0, 0.182223107579),\n",
    "            (0.,0.,0., 0,1,0, 0,1,0, 0.182223107579),\n",
    "            (0.,0.,0., 0,0,1, 0,0,1, 0.182223107579),\n",
    "\n",
    "            (0.,0.,0., 2,0,0, 2,0,0,  0.223223306785),\n",
    "            (0.,0.,0., 0,2,0, 0,2,0,  0.223223306785),\n",
    "            (0.,0.,0., 0,0,2, 0,0,2,  0.223223306785),\n",
    "\n",
    "            (1.,2.,3., 1,0,0, 1,0,0, -5.63387712455e-06),\n",
    "            (1.,2.,3., 0,1,0, 0,1,0, -0.000116463120359),\n",
    "            (1.,2.,3., 0,0,1, 0,0,1, -0.000301178525749),\n",
    "\n",
    "            (1.,2.,3., 2,0,0, 2,0,0, 0.00022503308545040895),\n",
    "            (1.,2.,3., 0,2,0, 0,2,0, 0.0006102470883881907),\n",
    "            (1.,2.,3., 0,0,2, 0,0,2, 0.0013427831014563411),\n",
    "\n",
    "            (0.,0.,0., 1,1,0, 1,1,0, 0.0136667330685),\n",
    "            (0.,0.,0., 0,1,1, 0,1,1, 0.0136667330685),\n",
    "            (0.,0.,0., 1,0,1, 1,0,1, 0.0136667330685),\n",
    "\n",
    "            (3.,2.,1., 1,1,0, 1,1,0, 5.976771621486971e-5),\n",
    "            (3.,2.,1., 0,1,1, 0,1,1, 1.5742904443905067e-6),\n",
    "            (3.,2.,1., 1,0,1, 1,0,1, 4.00292848649699e-6)\n",
    "        ]\n",
    "\n",
    "        val1 = vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,M)\n",
    "        val2 = vrr(cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,\n",
    "            aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,M)\n",
    "        @assert isapprox(val1,val2)\n",
    "        @assert isapprox(val1,result)\n",
    "        val3 = vrr_iter(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,M)\n",
    "        val4 = vrr_iter(cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,\n",
    "            aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,M)\n",
    "    end\n",
    "    println(\"test_vrr successfully completed\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "help_vrr_iter (generic function with 1 method)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function help_vrr(ax,ay,az,na,aI,aJ,aK,aexpn,bx,by,bz,nb,bexpn,\n",
    "    cx,cy,cz,nc,cI,cJ,cK,cexpn,dx,dy,dz,nd,dexpn,M)\n",
    "    # Helper function for vrr, since I for some reason changed the call signature\n",
    "    # from the python version. This version has the same signature as the python version\n",
    "    # and calls the Julia version\n",
    "    na*nb*nc*nd*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "        cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,M)\n",
    "end\n",
    "\n",
    "function help_vrr_iter(ax,ay,az,na,aI,aJ,aK,aexpn,bx,by,bz,nb,bexpn,\n",
    "    cx,cy,cz,nc,cI,cJ,cK,cexpn,dx,dy,dz,nd,dexpn,M)\n",
    "    # Helper function for vrr, since I for some reason changed the call signature\n",
    "    # from the python version. This version has the same signature as the python version\n",
    "    # and calls the Julia version\n",
    "    na*nb*nc*nd*vrr_iter(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,\n",
    "        cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,M)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test_vrr successfully completed\n"
     ]
    }
   ],
   "source": [
    "test_vrr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Timing comparisons\n",
    "Python/C times are from the int2e-speed.ipynb notebook in the same folder.\n",
    "\n",
    "## TODO\n",
    "* Get Python values for integrals for comparison\n",
    "* Fix vrr_iter calls\n",
    "* Look for time improvements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## s orbitals on the same center:\n",
    "9 usec time Julia compared to 0.75 usec in C and 49 usec in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000013 seconds (51 allocations: 2.594 KiB)\n",
      "  0.000011 seconds (51 allocations: 2.594 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4.373354581904758"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#                xyza   na   lmna   aa   xyzb    nb ab   xyzc      nc  lmnc   ac  xyzd   nd ad  M\n",
    "@time help_vrr(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0)\n",
    "@time help_vrr(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000124 seconds (81 allocations: 3.531 KiB)\n",
      "  0.000121 seconds (81 allocations: 3.531 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4.373354581904758"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time help_vrr_iter(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0)\n",
    "@time help_vrr_iter(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## sp orbitals on the same center:\n",
    "0.95 usec C, 58 usec Python, 13 usec Julia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000024 seconds (137 allocations: 7.328 KiB)\n",
      "  0.000024 seconds (137 allocations: 7.328 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#           xyza   na  lmna   aa  xyzb   nb ab  xyzc   nc  lmnc   ac  xyzd   nd ad  M\n",
    "@time help_vrr(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 1,0,0, 1., 0.,0.,0., 1.,1., 0)\n",
    "@time help_vrr(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 1,0,0, 1., 0.,0.,0., 1.,1., 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000130 seconds (95 allocations: 3.781 KiB)\n",
      "  0.000139 seconds (95 allocations: 3.781 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time help_vrr_iter(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 1,0,0, 1., 0.,0.,0., 1.,1., 0)\n",
    "@time help_vrr_iter(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 1,0,0, 1., 0.,0.,0., 1.,1., 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## sp integrals, 4 different centers\n",
    "1.2 usec C, 66 usec Python, 23 Julia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000024 seconds (135 allocations: 7.281 KiB)\n",
      "  0.000023 seconds (135 allocations: 7.281 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.4746004027185156"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time help_vrr(0.,0.,0., 1., 0,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,0,0, 1., 1.,1.,0., 1.,1., 0)\n",
    "@time help_vrr(0.,0.,0., 1., 0,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,0,0, 1., 1.,1.,0., 1.,1., 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000124 seconds (93 allocations: 3.734 KiB)\n",
      "  0.000119 seconds (93 allocations: 3.734 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time help_vrr_iter(0.,0.,0., 1., 0,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,0,0, 1., 1.,1.,0., 1.,1., 0)\n",
    "@time help_vrr_iter(0.,0.,0., 1., 0,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,0,0, 1., 1.,1.,0., 1.,1., 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## pd integrals, 4 different centers\n",
    "3 usec C, 95 usec Python, 32 usec Julia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000039 seconds (747 allocations: 40.594 KiB)\n",
      "  0.000032 seconds (747 allocations: 40.594 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-0.01222864716176887"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time help_vrr(0.,0.,0., 1., 1,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,1,0, 1., 1.,1.,0., 1.,1., 0)\n",
    "@time help_vrr(0.,0.,0., 1., 1,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,1,0, 1., 1.,1.,0., 1.,1., 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
